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ABSTRACT 

We present a series of numerical simulations aimed at understanding the nature and origin of 
turbulence in coronal loops in the framework of the Parker model for coronal heating. A coronal loop 
is studied via reduced magnetohydrodynamics simulations in Cartesian geometry. A uniform and 
strong magnetic field threads the volume between the two photospheric planes, where a velocity field 
in the form of a ID shear fiow pattern is present. Initially the magnetic field that develops in the 
coronal loop is a simple map of the photospheric velocity field. This initial configuration is unstable 
to a multiple tearing instability that develops islands with X and O points in the plane orthogonal 
to the axial field. Once the nonlinear stage sets in the system evolution is characterized by a regime 
of MHD turbulence dominated by magnetic energy. A well developed power law in energy spectra is 
observed and the magnetic field never returns to the simple initial state mapping the photospheric 
fiow. The formation of X and O points in the planes orthogonal to the axial field allows the continued 
and repeated formation and dissipation of small scale current sheets where the plasma is heated. We 
conclude that the observed turbulent dynamics are not induced by the complexity of the pattern that 
the magnetic field-lines footpoints follow but they rather stem from the inherent nonlinear nature of 
the system. 

Key words: magnetohydrodynamics (MHD) — Sun: corona — Sun: magnetic topology — turbulence 
Online-only material: animations at http://www.df. unipi.it/^rappazzo/shear/| {non permanent link) 
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1. INTRODUCTION 



In recent papers fRapp azzo et al.ll2QQ7l . 12008^ we have 
described reduced magnetohydr odynamics (RMHD) sim- 
ulations of the Parker problem (|Parkerlll972l I1988L 11994 ) 
for coronal loops in Cartesian geometry. We have shown 
that the system develops small scales, organized in cur- 
rent sheets elongated in the direction of the DC mag- 
netic field, through an MHD turbulent cascade, and that 
a well defined power law spectrum is developed for total 
energy. The energy spectra develop steep slopes [seen 
also in the similar sim ulations of a line-tied RMHD sys- 
tem by lDmitruk et al.l (j2QQ3i) ] with spectral indices going 
from the classical —5/3 Kolmogorov spectrum up to al- 
most —3. The energy spectral indices (slopes) have a 
bearing on the heating rate of the boundary-forced sys- 
tem. The heating rates are significantly increased for 
realistic values of the magnetic field intensity and loop 
lengi:th compared to the scaling laws with fixed indices 
(jDmitruk fc G 6mez"1999). the last being recovered when 
a Kolmogorov-like spectrum develops (i.e. for weak fields 
or long loops). 

In the published simulations we have always used pho- 
tospheric velocity patterns made up of large spatial scale 
projected convection cells {large-scale eddies), mimick- 
ing disordered photospheric motions. As a consequence 
the magnetic field developing in the corona was not in 
equilibrium, but dynamical evolution occurred from the 
outset with a time-scale set by the interplay of forcing 
and nonlinear ity. In this paper, in order to clarify the 
origin of the MHD turbulent dynamics found in our pre- 



vious work, we explore the dynamics of coronal loops 
system when a shear velocity fiow stirs the footpoints 
of the magnetic field-lines. It is commonly thought that 
the topology of the photospheric driver should strongly 
infiuence the dynamics of a coronal loop, and that the 
magnetic field-lines anchored to the photospheric planes 
should passively follow their footpoints motions. In this 
picture the electric currents should develop along neigh- 
boring field-lines whose footpoints have a relative shear 
motion. In particular it might be argued that in our pre- 
vious simulations the turbulent dynamics of the coronal 
medium originated from the "complexity" of the photo- 
spheric forcing patterns, with their large-scale eddies and 
stagnation points. 

The simple unidimensional shear forcing velocity used 
here allows a very clear-cut numerical experiment: if the 
field-lines were to passively follow the imposed footpoint 
motion only a sheared magnetic field would develop in- 
side the volume. The dynamics o f these curre nt layers 
are subject to tearing instabilities (jFurth et al.| [T963). If 
this were the physical process at work the topology of 
the magnetic field would remain a mapping of the forc- 
ing velocity pattern, periodically disrupted by tearing- 
like instabilities when the shear grew beyond a threshold 
amount. Simulations show that initially the magnetic 
field is sheared and a tearing instability develops, but af- 
terwards turbulent dynamics similar to those with more 
complex boundary patterns develop, clarifying that tur- 
bulence does not stem from a direct infiuence of the pho- 
tospheric velocity pattern but it is due to the inherent 
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nonlinear properties of the system. 

The Parker Scenario for coronal heating has been 
the subj ect of int e nse re se arch. Both an al ytical 
ISturroc k^ Uchidal flMl Ballegooiienl Il986l: 

Antiocho^ 19871: 'Berger 1991; Gomez & Ferro-Fontan 
19921; ICowlev et al. 1997; Ng & Bhattachariee 1998; 
Uzdenskv 2007; Janse & Low 2009; Low & Janse 2009; 
Berger fc Asgari- Targhi J009; Alv & Amari 2010) and 



numerical (Mikic et al.l 11989: Longcope & Sudan 1994: 
Hendrix fc Van Hoven 1996; Galsgaar d fc Nordlundl 
iH; [Huang fc Zweibel 2009; Huang et al.l I2010D 
investigations have been carried out. Given the 
complexity of such boundary forced field-line tan- 
gling systems simplified models have been developed, 
including 2D incompressible MHD mode ls with ma^- 
netic forcing (^inaudi et al. 1996; Dmi truk fc Gome; 
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19971: IG eorgouhs et aL| ll998i : iDmitruk et al.l 
Einaudi fc Vel h 1999.) and shell models (i Nigro et al .l 
2004; Buchli n fc^ellil]2007[ ). Insights gleaned from these 
models include the important result that dissipation in 
forced MHD turbulence occurs in the form of "bursty" 
events with well defined power law distributions in 
energy release, peak dissipation, and duration, in a way 
reminiscent of, and consistent with, the distribution of 
fiares in the solar corona. 

An analytical model of a forced system very simi- 
lar to the simulati on presented here was proposed by 
iHevvaerts fc PriestI ([1992), and recently extended to the 
anisotropic turbulence regime by Bigot et al. (2008). 
They started with same MHD system threaded by a 
strong axial magnetic field in cartesian geometry and 
apply at the top and bottom boundaries two ID ve- 
locity fields of opposite direction and assumed that the 
sheared structure that develops in the corona then dis- 
sipates via an effective "turbulent resistivity" provided 
by a cascade, so that a dissipative equilibrium is set 
up in which shearing is balanced by slippage provided 
by the turbulence. This amounts essentially to a one - 
point closure model of MHD turbulence (Biska nipll2003f ). 
where turbulence acts only on very small-scales while 
the large-scales remain laminar and indeed with the 
same large-scale magnetic structure. They then use 
the eddy-dampe d quasi-normal Markovian approxima- 
tion (EDQNM) (jPouquet et al.l I1976D to estimate the 
effective cascade and dissipation for the given driving 
shear, and this allow them to develop a heating theory 
in which the only free parameter is the equivalent Kol- 
mogorov constant. In contrast, our simulations will show 
that nonlinearity cannot be neglected even at the large- 
scales, so that the turbulent self-consistent state has lit- 
tle resemblance to the imposed shear fiow. As a result 
[Heyvaerts fc Priest (1992) overestimate the heating rate 
as laminar dynamics would lead to a higher energy in- 
jection (Poyntin g; fiux). 

More recentlv iDahlburg et al! (|2005L l2009f ) have pro- 
posed the so-called "secondary instability" as a lead- 
ing mechanism operating in the Parker Scenario, re- 
sponsible for the rapid release of energy. In their 
view, disruption on ideal time-scale must arise after 
some time while slow quasi-steady reconnection allows 
magnetic energy to continue to accumulate in the sys- 
tem. In their view the system evolution may be de- 
scribed by a sequence of equilibria destabilized by mag- 
netic reconnection. The bulk of numerical simulations 



performe d by [E inaudi et al.l (1 19961); IDmitruk fc GomezJ 
(1997); IGeorg ouhs et al. (1998); Dmitruk fc Gomez 
(d999i): [Ein audi fc Velh (1999|); ,Dmitruk et al. (2003); 
iRappazzoe t al. (2007, 2008) has proven that the system 
does not evolve through a sequence of equilibria, rather 
more complex dynamics develop. 

The initial setup of the simulation presented in 
([Dahlburg et al.ll2QQ9f ) is also very similar to the one im- 
plemented in this paper but, besides the lower resolution 
and therefore the higher infiuence of numerical diffusion, 
the time interval for which they advance the equations is 
too short compared with coronal loops and active region 
time-scales. These leads them to claim as representa- 
tive of the dynamics what is actually a transient event 
taking place only during the early stage of the dynam- 
ics, in our case the multiple tearing mode current sheet 
collapse. This evolution is not generic, but is only rep- 
resentative of a small class of very symmetric boundary 
velocity patterns, those which admit coronal equilibria 
at all times. We will return to this question and a more 
detailed discussion in the conclusion. 

The paper is organized as follows. In § [2] we describe 
the basic governing equations and boundary conditions, 
as well as the numerical code used to integrate them. In 
§ [3] we discuss the initial conditions for our simulations 
and briefiy summari ze the linear stage dyna mics more ex- 
tensively detailed in 'Rappazzo et al.' ('20 08[). while in_§j4] 
we outline the main points of Heyvaerts fc PriestI (jl992[ r 
relevant to this work. The results of our numerical sim- 
ulations are presented in § [Sj while the final section is 
devoted to our conclusions and discussion of the impact 
of this work on coronal physics. 

2. GOVERNING EQUATIONS AND BOUNDARY 
CONDITIONS 

We model a coronal loop as an axially elongated Carte- 
sian box with an orthogonal cross section of size i and 
an axial length L embedded in an homogeneous and uni- 
form axial magnetic field Bq = aligned along the 
z-direction. Any curvature effect is neglected. 

The top and bottom plates {z = and L) represent 
the photospheric surfaces where we impose, as bound- 
ary conditions, velocity patterns mimicking photospheric 
motions. Along the x and y directions periodic boundary 
conditions are implemented. 

At the top plate z = L we impose a sinusoidal shear 
flow with wavenumber 4 

27r 



u (x, y) = sin I 4 -— x + 1 



(1) 



At the bottom plate z = we generally impose a vanish- 
ing velocity 

u«(x,^)=0, (2) 

except in one simulation [run F (table [1])] where, in or- 
der to compare with previous simulations implement- 
i n^ a vortical velocity forc ing applied at both plates 
([Rappazzo et al.|[2007ll2008[ ). the reversed pattern of the 
top plate ([!]) is applied 



(x, y) = — sin ^4 ^ X + 1 



(3) 



Here Gy is the unitary vector directed along the y direc- 
tion, while the flow is sheared along x. 
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The dynamics are integrated, a s in our previous works, 
using the equ ations of RMHD (iKadomtsev fc Pogutsd 
[1971 IStraussI [1975 [Montgomery' 198 2"), which are weh 
suited for a plasma embedded in a strong axial magnetic 
field. In dimensionless form they are given by: 



the scale A associated with the dissipative terms used in 
equations ([I])-© is given by 



(b^ • V^)b^ +CA 



. (-l)-vr'„,, (4) 



dz 



+ (u^-V^)b^=(b^.Vju^+CA 



dt 



• = 0, • b^ = 0, 



dz 



(-1) 



n+l 



RCn 



(6) 



where and b^ are the velocity and magnetic fields 
components orthogonal to the axial field, p is the kinetic 
pressure. The gradient operator has components only in 
the perpendicular x-y planes 

while the linear term cx dz couples the planes along 
the axial direction through a wave-like propagation at 
the Alfven speed ca- Incompressibility in RMHD equa- 
tions f ollows from t he large value of the axial magnetic 
fields (lStrausslll976[) and they remain valid also for low 
P sys tems (|Zank fc Matthaeus|[T992[ : lBhattachariee et al.l 
1 19981 ) such as the corona. 

To render the equations nondimensional, we have first 
expressed the magnetic field as an Alfven velocity \b 
b/ v^47rpo], where po is the density supposed homogeneous 
and constant, and then all velocities have been normal- 
ized to the velocity ix* = 1 kms~^^ the order of magni- 
tude of photospheric convective motions. 

Lengths and times are expressed in units of the per- 
pendicular length of the computational box t = ^ and 
its related crossing time = As a result, the 

linear terms cx dz are multiplied by the dimensionless 
Alfven velocity ca = va/u*^ where va = Bq/^/A^pq is 
the Alfven velocity associated with the axial magnetic 
field. 

The majority of the simulations performed [specifically 
runs A-E (see table [T])] use a standard simplified diffu- 
sion model, in which both the magnetic resistivity 77 and 
viscosity v are constant and uniform. The kinetic and 
magnetic Reynolds numbers are then given by: 

Re = - , Re^ = 8 

where c is the speed of light, and numerically they are 
given the same value Re = Re^. In equations ([l])-(j5j) 
this case is realized for n = 1 with Re^ = Re. 

The index n is called dissipativity and forn > 1 the dis- 
sipative terms in ([I])-© correspond to so-called hyper- 
diffusion (Biskamp 2003). We use hyperdiffusion, with 
n = 4, only in runs F and G (table [1]) dedicated to study 
the energy spectra. Hyperdiffusion is used because, even 
with a grid of512x512 points in the x-y plane (the high- 
est resolution grid we used for the plane), the timescales 
associated with ordinary diffusion are small enough to 
affect the large-scale dynamics and render difficult the 
resolution of an inert ial range. The diffusive time at 



Re^ A'^. 



(9) 



Vf (5) 



For n = 1 the diffusive time decreases relatively slowly 
toward smaller scales, while for n = 4 it decreases far 
more rapidly. As a result for n = 4 we have longer 
diffusive timescales at large spatial scales and diffusive 
timescales similar to the case with n = 1 at the resolu- 
tion scale. Numerically we require the diffusion time at 
the resolution scale Xmin = 1/^, where N is the number 
of grid points, to be of the same order of magnitude for 
both normal and hyperdiffusion, i.e., 

^ Re^^Re,N'^^-'\ (10) 



Then for a numerical grid with N = 512 points that 
requires a Reynolds number Re^ = 800 with ordinary 
diffusion we can implement Re^ ~ 10^^ (table [1]), remov- 
ing diffusive effects at the large scales and allowing (if 
present) the resolution of an inertial range. 

We solve numerically equations (|1])-(|6|) written in terms 
of the potentia ls of the ort hogonal velocity and mag- 
netic fields [see 'Ra ppazzo et al.l (j2007|, 12008) for a more 
detailed description of the numerical code] in Fourier 
space, i.e. we advance the Fourier components in the x- 
and ^/-directions of the scalar potentials. Along the z- 
direction, no Fourier transform is performed so that we 
can impose non-periodic boundary conditions (§[3|), and 
a central second-order finite-difference scheme is used. In 
the x-y plane, a Fourier pseudospectral method is imple- 
mented. Time is discretized with a third-order Runge- 
Kutta method. 

3. INITIAL CONDITIONS AND LINEAR STAGE 

At time t = we start our simulations with a uniform 
and homogeneous magnetic field along the axial direction 
B = BoGz. The orthogonal component of the velocity 
and magnetic fields are zero inside our computational 
box = =0, while at the top and bottom planes a 
large-scale velocity pattern is imposed [([I])-([2j) or ([I])- ([3])] 
and kept constant in time. 

We briefiy summarize and extend to the shear forcing 
considered in this pa per the linear s ta^e analysis cov- 
ered in more detail inl Rappazzo et all (|2Q08). In general 
for an initial interval of time smaller than the nonlinear 
timescale t < Tnu nonlinear terms in equations ([4])- ([6]) 
can be neglected and the equations linearized. For sim- 
plicity we will at first neglect also the diffusive terms and 
consider their effect in the second part of this section. 
The solution during the linear stage for generic bound- 
ary velocity forcings, and respectively at the top 
and bottom planes z = L and 0, is given by: 

h^{x,y,z,t)= [u^(x,^)-uO(x,^)] — , (11) 

TA 

u^{x,y,z,t) = M^{x,y) | + uO(x,^) (l - |) , (12) 

where ta = L/va ''^^ the Alfven crossing time along the 
axial direction z. The magnetic field grows linearly in 
time, while the velocity field is stationary and the order 
of magnitude of its rms is determined by the boundary 
velocity profiles. Both are linear combinations {mapping) 
of the boundary velocity fields. 
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Considering the boundary conditions ([I])-© imposed 
in most of the fohowing simulations we obtain 

(x, y, z,t) = ^ sin ^4 ^ x + 1^ Gy , (13) 

u^{x,y,z,t) = - sin ( 4 — x + 1 1 Gy. (14) 

Both fields are a clear mapping of the shear velocity at 
the boundary, with the magnetic field increasing its mag- 
nitude linearly in time. 

When we shear the field-lines from both photospheric 
plates (2: = and L) using the forcing ([I])-(|3]) we obtain 
a very similar result, still a mapping but with different 
amplitudes: 

(x, y,z,t) = 2^ sin ^4 y x + 1^ Gy , (15) 
(x, y, z, t) = (2^ - 1) sin (^4 y X + 1^ Gy. (16) 

For a generic forcing the solution (pT]) - (p!2]) is valid only 
during the linear stage, while for t > Tni when the fields 
are big enough the nonlinear terms cannot be neglected. 

Nevertheless there is a singular subset of velocity forc- 
ing patterns for which the generated coronal fields (pT]) - 
(p!2]) have a vanishing Lorentz force and the nonlinear 
terms vanish exactly. This subset of patterns is charac- 
terize by having the vorticit y constant along the stream- 
lines (R appazzo et al.ll2QQ8f ). In this case the solutions 
are an exact solution at all times^ not an ap- 
proximation valid only for t < Tni- 

As can be proved by direct substitution the sheared 
forcing ([I])-© [and also ([I])-(|3])] is one of these degenerate 
patterns, and the solution (fl3|) - (p!4|) [or (p!5|) - (fT6|) ] is exact 
at all times. 

So far we have neglected the diffusive terms in the 
RMHD equations(j4j)-(j6]). Magnetic reconnection devel- 
ops even for small values of the resistivity, but in this 
section we are interested in the diffusive effects on the 
linear dynamics, i.e. when nonlinear terms are negligible 
or artificially suppressed (we will discuss such a case in 
our Conclusions). We now consider the effect of standard 
diffusion [case n = 1 in eqs. (j4])-(j5])] on the solutions (pT|) - 
(p!6|) : these are the solutions of the linearized equations 
obtained from (|4])-(|5]) retaining also the diffusive terms. 

In the linear regime^ as the magnetic field grows in 
time [([II]), ([I3]), ([ISj)], the diffusive term [V^ b^ cxb^/^^] 
becomes increasingly bigger until diffusion balances the 
magnetic field growth, and the system reaches a satu- 
rated equilibrium state. Considering diffusion, the mag- 
netic field will evolve as 



rate J will then be given by 



h^{x,y,z,t) = [\i^{x,y) - M^{x,y) 

TR 



TA 



1 



exp 

TR 



(17) 



The diffusive timescale associated with the Reynolds 
number Re is tr = ^^i?e/(27r)^ where ^c is the length- 
scale of the forcing pattern, that for the pattern ([I])-© 
is given by £c = where ^ is the orthogonal computa- 
tional box length. 
The total magnetic energy Em and ohmic dissipation 



E 



M 



2 Jv 



rrisat 



1 — exp 



1 — exp 

TR 



, (18) 
(19) 



where E^^ and J^^^ are the saturations value reached 



for t > 2ri^, whose values are given by: 



2L(87r)^ 



i(87r)2 • 



(20) 



Magnetic energy saturates to a value proportional to the 
square of both the Reynolds number and the Alfven ve- 
locity, while the heating rate saturates to a value that 
is proportional to the Reynolds number and the square 
of the axial Alfven velocity. In general the equations 
of RMHD are valid as far as the orthogonal magnetic 
field b^ is small compared to the dominant axial field 
Bo = BoBz. In particular incompressibility holds as 
far as the perturbed magnetic pressure can be neglected 
compared to that of the strong field b^ ^ ^a- There- 
fore the solutions found in this section are valid as far 
as the saturated values of the magnetic field satisfy the 
previous condition. In all the simulations presented here 
this condition is satisfied. 

4. EFFECTIVE DIFFUSIVITY: ONE-POINT 
CLOSURE MODELS 

Given the complexity of the Parker problem, simplified 
models have been derived. iHevvaerts fc PriestI ()1992l ) 
have developed an effective diffusivity model that is in 
effect a one-point closure model. In order to discuss the 
impact of our work on their results, we briefiy summa rize 
the relevant one-point closure theory ()Biskampll2QQ3f ). 

In order to investigate basic properties of the Parker 
model and compare them with observational constraints, 
such as the global heating rate and required energy fiux, 
the detailed dynamics of turbulent fluctuations (that are 
essential to determine how the individual field-lines are 
heated and hence how radiation is emitted) contain more 
informations than actually needed. 

It may then be attempted to split the velocity and 
magnetic fields into average and fluctuating parts: 



(B) 



u. 



(21) 



Incompressible MHD equations give for mean flelds the 
following equations: 



St(u) + (u)-V(u) 



-V(P) + (B)^V(B) 

-V- (uu-bb) + i/V^(u), (22) 



5t(B) + (u)-V(B) = (B)-V(u) 

+V X (S X b) +7/V^(B), (23) 

V • (u) = V • S = V • (B) = V • b = 0, (24) 

where symbols have the usual meaning, and in particular 
V and 77 are the microscopic viscosity and resistivity of 
the plasma. If it is possible to model the terms that 
contain the small-scale fluctuations then eqs. (|22]) - (|24|) 
allow to advance the mean flelds. 
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Table 1 

Summary of the simulations 



Run 



CA Ux X Uy X Hz 



forcing 



Re or Re4 



A 


200 


512 


X 


512 X 


200 


shear: 


t 


1 


800 


B 


200 


256 


X 


256 X 


100 


shear: 


t 


1 


400 


C 


200 


128 


X 


128 X 


50 


shear: 


t 


1 


200 


D 


200 


128 


X 


128 X 


50 


shear: 


t 


1 


100 


E 


200 


128 


X 


128 X 


50 


shear: 


t 


1 


10 


F 


200 


512 


X 


512 X 


200 


shear: 


t,b 


4 


1019 


G 


200 


512 


X 


512 X 


200 


vortex 


t,b 


4 


1019 



Note. — CA is the axial Alfven velocity, Ux x Uy x Uz is the 
numerical grid resolution, b or t in the forcing column indicates 
whether the forcing is applied at the bottom and/or top plates, 
n is the dissipativity, n = 1 indicates normal diffusion, n = 4 
hyperdiffusion. The next column indicates respectively the value 
oif the Reynolds number Re or of the hyperdiffusion coefficient Re4. 



This is a two-scale approach, and the average and 
fluctuating parts can be represented as the large-scale 
and small-scale flelds or, introducing a suitable cut-off 
wavenumber K (for instance the grid resolution if this is 
applied to numerical simulations), as the low- and high- 
pass filtered "lesser" and "greater" functions 



(B)=B< = ^Bke 



ikx 



k<K 



b = B> = 



ke 



ikx 



(25) 



(26) 



k>K 



The relevant quantities to be modeled in eqs. 
are the "turbulent" stress tensors 



Ri 



Sin 



-{uibj 



Ujbi). 



(27) 



Phenomenolog!:ical mo deling and numerical simulations 
(|Yoshizawalll99QL [1991 ) have shown that these stress ten- 
sor can be approximated with 



Ri 



Sin 



Vt 



(28) 
(29) 



where the coefficients Ut and rjt are called turbulent vis- 
cosity and resistivity, and for a plasma in coronal condi- 
tions have much higher values than u and rj. Inserting 
([Ml)-® in the equations for the mean fields (f22)) - (|23)) . 
we obtain a set of equations for the mean fields that has 
the same structure of incompressible MHD equations ex- 
cept that u ^ Uf and r] ^ rjt. 

Then for a system embedded in a strong axial field 
([Montgomery 1982) equations ([4])- ([6]) that we use to 
model the system are also a one-point closure model if the 
d issip ative coefficients are co nsidered as effective ones. 

iHevvaerts fc Priest! ()1992l ) use the results of an 
eddy damped quasi- normal Markovian approximation 
(iPouquet et al.| [l976[ ) to express the effective diffusivity 
coefficients Vt and rjt as a function of the energy flux e 
flowing along the inert ial range. 

Additionally they suppose that the ID boundary forc- 
ing velocity leads the large-scales to evolve in a laminar 
(flelds are ID and directed along the same direction of 
the forcing) steady state {dt = 0). Equations ([4])- ([6]) and 



their boundary forcing [of which our forcing ([I])- ([2]) is 
representative] are used to compute the flux of energy S 
entering the system due to the dragging of the fleld-lines 
footpoints. 

/S is a function of the system parameters and the ef- 
fective diffusivity coefficients Ut and r]t that are the only 
unknown variables in both energy ffuxes S and e. The 
solution of the problem is achieved requiring that the ffux 
S entering the system at the large scales is equal to the 
ffux e ff owing from the large to the small scales. 

As a matter of fact their laminar steady state coincides 
with the saturation linear state that we have computed 
in the previous section (§[3|). This was obtained neglect- 
ing the nonlinear terms in eqs. ([4])- ([5]) and retaining the 
diffusive terms. They obtain it in a similar way, by sup- 
posing that the induced flelds will retain the ID sym- 
metry of the forcing. Thus the nonlinear terms vanish, 
as can be proved by direct substitutions of the generic 
flelds = f{x) Gy, = g{x) Gy, with / and g generic 
functions. 

Our simulations investigate the large-scale dynamics 
and as we show in § [5] the large-scale flow is not laminar 
and it is steady only in a statistical sense. Turbulence 
cannot be conflned only to the small scales. We dis- 
cuss the implications for t he flndings and scaling laws of 
IHevvaerts fc PriestI (|1992r ) in section § [6] devoted to our 
conclusions. 

5. NUMERICAL SIMULATIONS 

In this section we present a series of numerical simu- 
lations, summarized in Table [TJ We will flrst describe 
the results of simulations A-E that model with different 
resolutions and associated different Reynolds numbers a 
coronal layer driven by the sheared velocity pattern ([1]) 
at the top plate z = L and a vanishing velocity ([2|) at the 
bottom plate z = 0. In all simulations the computational 
box has an aspect ratio of 10 with i = 1 and L = 10. 

5.1. Shear Forcing: run A 

We present here the results of run A, a simulation 
performed with a numerical grid of 512 x 512 x 200 
points, normal (n=l) diffusion with a Reynolds num- 
ber Re = 800. The Alfven velocity is va = 200 km 
corresponding to a ratio ca = VA/uph = 200. The total 
duration is 600 axial Alfven crossing times ta = L/va- 

We add to the system a perturbation (naturally present 
in the coronal environment) for the magnetic and velocity 
flelds inside the computational box at time t = 0. If there 
were no perturbation the system would follow the linear 
saturation curves plotted with dashed lines in Figures [T] 
and [21 as discussed in ^ 

We have not performed a full linear instability analysis, 
but we have used as perturba tions either an Orszag-Tang 
vortex (|Orszag Tand[T97l ). i.e. 



Su- 



Sh = e 



sin {27ry) Gx + sin (27rx) Gy 



sin {4:7ry) Gx + sin (27rx) Gy 



(30) 



(31) 



or a "white noise" (i.e. the value of the potentials associ- 
ated with the flelds are given in each grid point a random 
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100 200 300 t/t^400 500 600 



Figure 1. Run A: Magnetic (Em) and kinetic (Ek) energies as 
a function of time (r^ = L/va is the axial Alfven crossing time). 
The dashed curve shows the time evo lution of magnetic energy if 
the system were unperturbed [eq. (|18|) ]. 

number included between and 1), with different ampli- 
tudes e. We have computed the resulting growth rates 7 
with our nonlinear code. 

For both perturbations the value of the growth rates 
is high and similar to each other. The white noise has a 
slightly higher growth rate at ^ta ~ 0.69. This timescale 
is of the order of the Alfven crossing time r^, which 
corresponds to an ideal timescale. This implies that the 
Orszag-Tang vortex is close to the most unstable mode, 
which is selected from all the modes that are excited by 
the random perturbation. As we shall see in the following 
the instability is a multiple tearing mode. 

For each kind of perturbation, the smaller the ampli- 
tude the later the system becomes unstable. For the sim- 
ulation presented in this section we have used a "white 
noise" perturbation with an amplitude e = 10~^^, very 
small compared to the boundary imposed velocity eq. ([Tj) 
which is 

In Figures [l]l2] we plot the total magnetic and kinetic 
energies 

EM = \JdVhJ, EK = \JdVuJ, (32) 

and the total ohmic and viscous dissipation rates 

J=iJdVf, ^=^J<^V.\ (33) 

along with the saturation curves (p!8]) - (p!9|) . For smaller 
values of the perturbation amplitude e the system be- 
comes unstable sooner but the ensuing dynamics in the 
fully nonlinear stage are similar with same average values 
for all relevant physical quantities. 

As shown in Figures [l]l2] and [3] the shear velocity at 
the top boundary {z = 10) initially induces velocity and 
magnetic fields inside the volume that follow the linear 
behavior [eq. (p!2)) and (pT]) ]. The magnetic energy and 
the ohmic dissipation rate also follow initially the linear 
diffusive saturation curves eq. (p!8|) - (p!9|) . From eq. (|2Q]) 
with the values for this simulation we have that the mag- 
netic energy and ohmic dissipation would reach the sat- 
uration values 

E'^' = 1604, J"^' = 2533, (34) 

with a diffusion time tr ~ 25ryi, if the system were 
unperturbed. 

Figure [3] shows snapshots of the magnetic field-lines 
(orthogonal component b^^) and electric current at se- 
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Figure 2. Run A: Ohmic (J) and viscous (^2) dissipation rates 
as a function of time. The dashed curve is th e hnear saturation 
curve for the ohmic dissipation rate [eq. (|19|) ]. In the inset the 
integrated Foynting flux S is shown to dynamically balance the 
total dissipation. 



lected times in the mid-plane z = 5. Until time t ^ 7brA 
the velocity and magnetic fields in the volume are a 
mapping of the velocity at the boundary [eq. (p!2]) and 
(pTj)]. i.e. the sheared velocity at the boundary induces a 
sheared magnetic field in the volume. 

This magnetic configurati on is well known to be unsta- 
ble to tearing instabilities (Furt h et al.l ll963). and in fact 
around time t ~ 79rA a multiple tearing instability de- 
velops when the system transitions from the linear to the 
nonlinear stage. Magnetic field-lines reconnect in corre- 
spondence of X-points while the characteristic magnetic 
islands are formed. The randomness of the perturbations 
is refiected in the lack of symmetry of the reconnecting 
field. 

The magnetic energy (Figure [1]) accumulated up to this 
point is then released in a big burst of ohmic dissipation 
(see Figure [2]). Notice that correspondingly there is a 
peak in the viscous dissipation (i.e. the integrated vor- 
ticity), but it is smaller that the electric current peak 
as in the following dynamical evolution. This is in fact 
a magnetically dominated system also when considering 
only the magnetic fiuctuations with respect to veloc- 
ity fiuctuations , without taking into account the big 
axial field Bq. 

Up to this point the dynamics are not surprising, the 
magnetic field gets sheared and the shear increases in 
time. Such a configuration is very well known to be un- 
stable to reconnection, and then it is expected that a 
tearing- like instability should develop. 

What is surprising are the dynamics at later times. 
In fact as can be seen in Figures [I]l2] for t > 90 r a 
the system reaches a magnetically dominated statisti- 
cally steady state where integrated quantities fiuctuate 
around an average value. In particular velocity fiuc- 
tuations are smaller than magnetic fiuctuations, which 
in turn are small compared to the axial magnetic field 
((b^)V2/^Q - 0.027). Therefore the orthogonal mag- 
netic and velocity fields satisfy the RMHD ordering that 
requires them to be small compared to the axial field Bq. 

The energy power entering the system at the bound- 
aries as a result of the work done by photospheric motions 
on the footpoints of magnetic field-lines is given by the 
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Figure 3. Run A: Axial component of the current j (in color) and field-lines of the orthogonal magnetic field in the midplane (z = 5) at 
selected times covering the linear and nonlinear regimes up to t ~ 200 r^. During the linear stage (t ~ 60 r^) the orthogonal magnetic field 
is a mapping of the boundary shear velocity [eq. ([T])] . After the transition to the nonlinear stage due to a multiple tearing instability {t ~ 79 
and ^2ta) the topology of the magnetic field departs from the boundary velocity mapping and evolves dynamically in time {t ~ 90, 110 
and 200 ta). 



integrated Pointing flux 
S 



ca / da b , 

z=L 



(35) 



where is the photospheric forcing velocity ([T]). 

Dissipation rates and the Poynting flux also fluctuate 
around a mean value. In particular, as shown in the inset 
in Figure [21 the Poynting flux and the total dissipation 



(ohmic plus viscous) balance each other on the average, 
although on shorter timescales a lag between the two 
signals is present. In this steady state the energy that 
is injected into the system is then, on the average, com- 
pletely dissipated. 

The surprising feature is that the shearing in the mag- 
netic field inside the volume is not recreated, as shown 
in Figure [3l It is natural to think that after the first 
big dissipative event around t ~ 82ryi, the shear fore- 
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Figure 4. Side and top views of a snapshot of magnetic field-lines (top row) and current sheets (bottom row) at time r ~ 550 ta- The box 
has been rescaled for an improved visualization. Top: Field lines of the total magnetic field (orthogonal plus axial), and in the midplane 
(z = 5) field-lines of the orthogonal component of the magnetic field. In the side view we have superimposed in yellow the streamlines 
of the boundary forcing velocity in the top plate (z = 10), at the bottom we impose a vanishing velocity. Bottom: Two isosurfaces of 
the squared current j^. The isosurface at the value = 2.8 x 10^ is represented in partially transparent yellow, while red displays the 
isosurface with = 8 x 10^, well below the maximum value of the current at this time j^^x — ^-^ ^ . As is typical of current sheets, 
isosurfaces corresponding to higher values of p are nested inside those corresponding to lower values. The current sheets filling factor is 
small. 



ing velocity at the boundary would recreate over time 
a sheared magnetic field in the system that should then 
lead to another big dissipative event and so on. 

The dynamics of this system are in fact commonly 
approximated as a sequence of equilibria, each destabi- 
lized by magnetic reconnection. This approximation is 
attained by neglecting the velocity and kinetic pressure 
in the MHD equations whose solution is then bound to 
be a static force-free equilibrium. Our simulations con- 
firm that the system is magnetically dominated, and in 
particular magnetic energy is bigger than kinetic energy, 
as shown in Figure [T] where on the average Em ^ 61Ek. 
But the self-consistent evolution of the kinetic pressure 
and velocity, although small compared with the domi- 



nant axial magnetic field does not bind the system 
to force-free equilibria and allows the possible develop- 
ment of alternative dynamics. 

In the following sections we will analyze further as- 
pects of the dynamics and the spectral properties of the 
system. But first we illustrate the topology of the mag- 
netic and velocity fields, to understand why a sheared 
magnetic field is not recreated. 

5.1.1. Magnetic Field Topology and the Origin of 
Turbulence 

As shown in Figure [3] at time t ~ 79 reconnection 
starts to develop, enhancing the ohmic dissipation, that 
reaches a peak around t ^ 82 r a- This big dissipative 
event burns a large fraction of the magnetic energy pre- 
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Figure 5. Run A: Snapshots at time t = 350 of the axial component of the current j (in color) and field-lines of the orthogonal 
magnetic field in 3 distinct x-y planes: the plane z = 1 close to the bottom unforced boundary, the midplane z = 5, and in plane z = 9 
close to the forced boundary z = 10. The box lenght is L = 10. 



viously accumulated by the system (Figured]). But not 
all the magnetic energy is dissipated, which would oth- 
erwise bring the system to a configuration similar to the 
initial condition at t = 0. 

This dissipative event is in fact due to magnetic recon- 
nection, that during its evolution produces a component 
of the magnetic field along x, the cross-shear direction, 
forming magnetic islands [see Figure [3] at times t ^ 79 
and 82 ta] • In fact around t ~ 90 r^, at the end of the big 
dissipative event, the topology of the orthogonal compo- 
nent of the magnetic field is characterized by magnetic 
islands. Naturally the Lorentz force does not vanish now 
and the vorticity is not constant along the streamlines. 
As typical of magnetic reconnection vorticity forms mis- 
aligned quadrup olar structures around current sheets [see 
[Rappazzo et al.l (12008)]. 

Although the forcing velocity at the boundary is always 
a shear [eqs. ([I])-(j2j)], nonlinear terms do not vanish as 
they do during the linear stage for t < 79 r a- When 
they vanish magnetic energy can be stored without get- 
ting dissipated (see 33]). But now nonlinearity can redis- 
tribute along the cross-shear (x) direction part of the en- 
ergy associated with the shear-aligned (^/-oriented) field 
along which the forcing injects energy, and continuously 
cascades to lower scales as described in the following sec- 
tions. 

The three-dimensional structures are shown in Fig- 
ure |4j Although the magnetic energy dominates over 
the kinetic energy, the ratio of the rms of the orthogonal 
magnetic field over the axial dominant field Bq is quite 
small. For ca = 200 it is ~ 3%, so that the average in- 
clination of the magnetic field-lines with respect to the 
axial direction is just ~ 2°, it is only for lower value of 
Ca that th is ratio increases and the angle increases ac- 
cordingly (jRappazzo et al.ll2008f ). The field-lines of the 
total magnetic field at time 550 are shown in Figure |4] 
{top row). The computational box has been rescaled for 
an improved viewing, and to attain the original aspect 
ratio, the box should be stretched 10 times along the ax- 
ial direction. The magnetic topology for the total field is 
quite simple, as the lines appear slightly bent. 

Figure^ {bottom row) also shows a view from the side 
and the top of the 3D current sheets at time 550 r^. The 
current sheets, elongated along the axial direction, look 
space filling when watched from the side of the compu- 
tational box, but the view from the top shows that the 
filling factor is actually small, as they are almost 2D 



structures. 

So far we have analyzed the topology of the field-lines 
only in the mid-plane z = 5. In Figure [5] we show the 
current density and magnetic field lines of in the mid- 
plane {z = 5) and in two other x-y planes close to the 
boundaries {z = 1 and 9, the axial length L = 10). The 
behavior is similar at different heights although in the 
plane z = 9, closer to the forced boundary z = 10, the 
topology of the field appears to be affected to some extent 
by the sheared velocity forcing ([T]) directed along the 
y direction. The field-lines in fact show a small alignment 
directed along y close to the boundary z = 10. 

The influence of the boundary forcing over the mag- 
netic field can be expressed quantitatively through the 
correlation between the magnetic field in the plane z 
and the boundary forcing velocity [eq.([l])]: 



Cor [b,,u^] {z) = 



JJdxdy b^ • 



JJ dxdy b^ . jj dxdy (u^)^ 



(36) 



In Figure [6] we plot the correlation as a function of the 
axial coordinate z at selected times. In the linear stage, 
until time slightly bigger than t = TGta whereafter the 
system transitions to the nonlinear stage (Figures [H [21 
and [3]), the magnetic field is a mapping of the boundary 
velocity therefore as expected the correlation is equal to 
1. In fact for the simulation presented in this section 
(run A), for which we have imposed the shear velocity 
profile ^ at the top plate z = 10 and a vanishing velocity 
at the bottom plate z = 0, the magnetic field in the 
linear stage is given by eq. (p!3|) [or (p!7|) with = 
including diffusion] therefore the correlation is 1 as b^ is 
proportional to the boundary velocity u^. 

Next as the system transitions to the nonlinear stage 
releasing most of the accumulated magnetic energy the 
correlation between the magnetic field and the boundary 
forcing velocity decreases swiftly, at a faster pace the 
farther from the forced boundary z = 10, as shown by 
the curves at times 82.19 < t/rA < 88.28. 

The correlation during the fully nonlinear stage is 
shown with color lines at 10 selected times separated 
by At = 40rA in the interval 200 < t < 600ta. 
The correlation vanishes near the bottom boundary and 
then grows almost linearly with z up to ~ 0.6 at the top 
boundary. As expected the correlation is bigger near the 
forced boundary and fades towards the interior of the 
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Figure 6. Correlation [eq. p6|) ] between the magnetic field 
and the boundary forcing velocity [applied at the boundary 
z = 10, see as a function of the axial coordinate z at se- 
lected times. Correlation is equal to 1 during the linear stage 
{t < 76 r a) and decreases as the system transitions to the non- 
linear stage (82.19 ^ t < 88.28 ry^). The different colors show 
the correlation during the fully nonlinear stage at 10 different times 
separated by At = 40 ta in the interval 200 ta ^ t < 600 ta • 

computational box. The magnetic field is overall weakly 
correlated with the forcing velocity, and at most reaches 
a mild correlation close to the forcing boundary, but it 
is never close to a strong correlation Cor = 1. 

In Figure [71 we show the 2D averages in the x-y planes 
of the magnetic and velocity fields and of the ohmic 
dissipation / R plotted as a function of z at differ- 
ent times. The behaviour is very similar to our previ- 
ous simulations with different {vortical) forcing patterns 
([Rappazzo et al. 2008). These macroscopic quantities 
are smooth and present almost no variation along the 
axial direction. The velocity must approach its bound- 
ary values at z = and 10, and in the volume grows 
to values higher than the boundary velocity ul- In fact 
also the velocity inside the volume is not a mapping of the 
boundary forcing but develops self-consistently, in partic- 
ular the plasma jets at reconnection locations contribute 
too to its average. The reason of the overall smooth be- 
havior of these quantities is that every disturbance or 
gradient along the axial direction is smoothed out by the 
fast propagation of Alfven waves along this direction; 
their propagation time ta is in fact the fastest timescale 
present (in particular faster than the nonlinear timescale 
< ^ni)i and then the system tends to be homogeneous 
along this direction. 

5.2. Dissipation versus Reynolds Number and 
Transition to Turbulence 

Simulations A-E differ only for the value of the 
Reynolds number (and corresponding grid resolution), 
all other parameters are the same including the ampli- 
tude of the perturbations. In Figure [8] we plot the total 
dissipation for the 5 simulations as a function of time. 

The simulation described in the previous paragraph 
had Re = 800, after the first big dissipative peak its 
signal displays a complex temporal structure that arises 
from the underlying turbulent dynamics. Decreasing the 
value of the Reynolds number to Re = 400 the structure 
of the signal has a simpler structure, with an almost si- 
nusoidal form of different amplitude and period ~ Sta in 
some time intervals. 

The dotted lines represent the linear diffusive behavior 
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Figure 7. Two-dimensional averages in the x-y planes of the 
ohmic dissipation j^/Re, the magnetic field b^, and the veloc- 
ity field , as a function of the axial coordinate z. The different 
colors represent 10 different times separated by At = AOta in the 
interval 200 ta < t < GOOta- 

described in § [3l These solutions [eq. (fT9|) ] are obtained 
when the nonlinear terms can be neglected. This can 
happen either because they vanish exactly due to the 
symmetry of the forcing velocity (as discussed in § [3]) or 
equivalently when nonlinear terms are depleted by diffu- 
sion^ i.e. when diffusion dominates the dynamics. 

At Re = 200 diffusion affects the system substantially, 
the instability takes longer to develop, and afterward 
the actual signal and linear saturation curve differ lit- 
tle. There is no first big dissipative peak, for a long time 
up to t ~ 550 the signal follows the linear saturation 
curve (p!9|) . Afterwards it departs slightly displaying a 
sinusoidal behavior of very small amplitude. 

At lower Reynolds number {Re = 100 and 10 are 
shown) diffusion dominates the dynamics, and no small 
scales is formed, while ohmic dissipation and energy fol- 
low the curve (p!8|) - (p!9|) describing the diffusive equilib- 
rium that is formed, as nonlinear terms are completely 
depleted. 

It is very interesting to notice that total dissipation 
for simulations A, B and C with Re = 800, 400 and 
200 substantially overlap each other. As typical with 
spectral numerical codes, which use Fourier transforms 
to compute derivatives, dissipation due to numerical im- 
plicit diffusion is very small. We have checked that the 
energy conservation equation, which can be derived from 
eqs. 0])-(|6]), is numerically verified within a very small 
error: it is around 2% for the lower resolution simula- 
tions and decre ases below 0.5% at higher resolutions (see 
lRappazzo"2006', Fi gure 5.5). Therefore in the simulations 
presented here the dissipation rates are substantially due 
to the explicit dissipative terms present in equations (|4])- 
([5]), while the diffusion due to the numerical schemes is 
negligible. 

We had already seen this b ehavior in our previous sim- 
ulations with vortex forcing (jRappazzo et al.ll2008l ). but 
in that case the boundary velocity was slightly differ- 
ent in each simulation. In fact it was built by a linear 
combination of Fourier modes with random amplitudes 
normalized to have the same rms, but the spatial pattern 
was different for every simulation as determined by the 
random amplitudes. In the simulations presented here 
the forcing is exactly the same for each simulation, as it 
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Figure 8. Total (Ohmic plus viscous) dissipation rates as a func- 
tion of time for different Reynolds numbers. The dashed curve s ar e 
the linear saturation curves for the ohmic dissipation rates [eq. (|19|) ] 
for the different cases. The overlap for Re > 200 suggests that total 
dissipation is independent from Reynolds number at higher values. 



is simply a single Fourier mode [eq. ([T])]. 

As the forcing is the same for each simulation the over- 
lap is more evident, and makes stronger the claim that 
total dissipation is independent of the Reynolds num- 
ber beyond a threshold. W e conjectured this hypothe- 
sis appazz o et al.l l2QQ8h because the energy flux en- 
tering the system (Poynting flux) and the energy trans- 
port from large to small scales due to the development 
of a turbulent dynamics are both independent from the 
Reynolds number, where the turbulent transport exhibits 
this property only for a sufficiently high value of the 
Reynolds number. 

This unfortunately does not imply at all that the ther- 
modynamical and radiative outcome is independent of 
the Reynolds number. In fact how fleld- lines are heated 
strongly depends on the dynamics and properties around 
the single current sheet. These become thinner and thin- 
ner at higher Reynolds numbers and the overall dynamics 
more chaotic. Hence the thermodynamical properties of 
the fleld-Hnes that cross the current sheets (elongated 
along the axial direction) and get so heated impulsively 
are all to be explored. 

An open question is whether the dissipation is inde- 
pendent of the Reynolds number also in the single cur- 
rent sheets, or their number and properties change to 
attain independence for the total dissipation. Research 
in this area is active (iLoureiro et a l. 2007; Lapenta 2008; 
Servidio et al."2009': 'Samtanev et al. 2009; Cassak et al. 
2009; Bhattachariee et al. 2009; Loureiro et al. 2009; 
Huang & Bhattachariee 2010; Servidio et al.i i2010„) and 
has already shown that at relatively high Reynolds num- 
bers reconnection departs the classic Sweet-Parker scal- 
ings in the MHD regime. Furthermore for a plasma in 
coronal conditions kinetic effects cannot be excluded a 
priori, in fact they might play an important role in dis- 
sipating energy through particle acceleration. Also the 
high value of the magnetic Prandtl number could have a 
bearing (Schekochihin et al. 2004). Nonetheless the total 
dissipation, integrated over the whole volume, is likely to 
not depend on the detailed small-scale dissipative mech- 
anism. In fact the energy injection rate [eq. (|37|) ] de- 
pends only on the large scale flelds and the transfer of 
this energy toward the small scales appears to be lo- 
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Figure 9. Transition to turbulence, total ohmic and viscous dissi- 
pation as a function of time form simulations A, B and C (displayed 
on the same scale for an equal interval of time). All the simula- 
tions implement ca = 200, but different Reynolds numbers, from 
Re = 200 up to 800. For Reynolds numbers lower than 100, the sig- 
nal is c ompletely flat following exactly the linear saturation curve 
eq. (|19|) . At higher Reynolds numbers, smaller temporal structures 
are present displaying a transition to turbulence. 



cal (Rappa zzo fc Vellill2QlQf ), i.e. it is determined by the 
flelds at neighboring large scales, thus making also the 
transfer energy rate toward the small scales independent 
of the Reynolds number (provided it is beyond a minimal 
threshold to have scale separation). 

In Figure [9] we show a close-up of total dissipation for 
an equal interval of time At = 600 ta on the same scale 
in order to highlight their temporal structures. At higher 
Reynolds number smaller time frequencies are p resent, 
a clear indication of a transition to turbulence (|Frischl 
1995). For Reynolds numbers lower than 100, the signal 
is completely flat following exactly the linear saturation 
curve eq. ^T^ . 

5.3. Spectral properties 

In order to study the spectral properties of the system 
and compare them with those of previous simulations we 
have performed a new simulation where the sheared forc- 
ing is applied on both the top and bottom plates (run F) 
with reversed direction [eqs. ^ and (|3|)]. We compare 
these results with those obtained from a previous simu- 
lation (run G), that has all the same parameters, except 
that on the two boundary planes a large-scale "vortex- 
like" velocity pattern with the same rms ((u^) = 1/2) is 
applied. Both simulations have been performed with a 
numerical resolution of 512 x 512 x 200 grid points, and 
hyperdiffusion with dissipativaty n = 4 and R4 = 10^^. 

In Figure [To] we plot the energy spectra obtained from 
runs F and G. They are very similar to each other. Both 
the velocity and magnetic flelds develop inertial ranges 
following power laws, and overlap each other. For both 
runs the spectral index of the kinetic spectrum (~ —0.5) 
is much smaller than that of the magnetic energy (~ 
—2.1), that is steeper than kolmogorov (—5/3). Also 
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Figure 10. Kinetic (Ek) and Magnetic (Em) Energy spectra 
as a function of the orthogonal wavenumber for simulations F 
and G. The spectra are very similar whether a vortex-like or shear 
velocity pattern stirs the footpoints of the magnetic field-lines. 

the kinetic energy has a lower value than the magnetic 
energy, as already noticed for the integrated quantities 
(Figured]). 

The energy that is injected into the system for unit 
time is the integrated Poynting flux 

S = CA [ da(u^ • ) - CA / da(u^ • ) (37) 

Jz=L Jz^O 

where and are the imposed velocity patterns at 
the top and bottom planes. 

For run G, we excite all wavenumbers 3 < < 4, 
while for run F we excite only one Fourier component as 
= — = sin (Sttx + 1) Gy, i.e. we are injecting en- 
ergy in the system only at iiin = 4 • 27r Gx, the wavenum- 
ber 4 along X. This can be noticed also in Figure [TOl 
where the kinetic spectrum for run G at n = 3 is higher 
than for run F, as part of the energy is injected also at 
n = 3 in the vortical case. 

The lower level for the kinetic spectrum is due to the 
boundary conditions that roughly set the value or the 
velocity at the injection wavenumbers inside the volume. 
In the simple linear case this is given by eq. (fT2]) . In the 
shear case we would have Ek{^) = 1/2 • V • (u^) =2.5 
(V = 10 is the volume) in the linear regime, and from 
Figure [10] we notice that also in the nonlinear regime 
Ek{^) ~ 2.5. On the other hand the magnetic field grows 
linearly in time [eq. (pT]) ] until a balance is reached be- 
tween the energy fiux that is injected at this scale and 
the fiux of energy fiowing towards smaller scales through 
a turbulent cascade. 

The magnetic energy spectra of the two simulations are 
slightly different at the large scales with < 5. The 
large scale dynamics is in fact slightly different in the two 
cases. 

In the vortical case (run G) energy is injected in all 
modes with wavenumbers 3 < < 4 that then cas- 
cades toward smaller scales. In the shear case (run F) 
energy is injected only at o ne wav enumber = (4,0). 
We have already noticed in ? l5.1.1l that although we con- 
tinue shearing the footpoints of the field-lines with our 
ID forcing [eq. ([T])] in the nonlinear stage the orthogonal 
magnetic field is organized in magnetic islands (Figure[3]), 
so that it is no longer a mapping of the boundary veloc- 
ity. Although energy is injected only in the wavenumber 
4 along X, energy is then redistributed by the nonlinear 



terms also to modes with wavenumbers along y at the 
large scales, and a small inverse cascade is present as in 
run G. This is the basic mechanism by which magnetic 
islands are sustained throughout the simulation in the 
nonlinear stage. 

6. CONCLUSIONS AND DISCUSSION 

In this paper we have investigated the dynamics of the 
Parker problem for the heating of coronal loops when 
the footpoints of the magnetic field-lines are stirred by a 
ID shear velocity pattern at the photosphere- mimicking 
boundary, and compared these results with those previ- 
ously obtained when a mo re complex ^Vortex- lik e" ve- 
locity pattern was imposed (iRappazzo et al.|[2008[ ). This 
very simple forcing is ideal to investigate the origin of tur- 
bulence in coronal loops and the influence of the bound- 
ary velocity forcing on the dynamics of the system. 

We will also co mpare our results with those of 
iHevvaerts^fc Pries t (|1992D and of the more recent sim- 
ulations of !.Dahlburg et al.l (|2005L [2009'). 

In summary, the main results presented in this paper 
are the following: 

1. Initially the sheared velocity forcing induces a sheared 
perpendicular magnetic field inside the volume. The 
resulting curren t layers are know n to be unstable to 
tearing modes (»Furth et al.l 119631 ). In fact when the 
system transitions from the linear to the nonlinear 
stage it is due to a multiple tearing instability, as 
shown in Figure [3] But once the system has become 
fully nonlinear the dynamics are fundamentally dif- 
ferent. As the nonlinear terms no longer vanish they 
now do transport energy from the large to the small 
scales where in correspondence of the X-points non- 
linear magnetic reconnect ion takes place, without go- 
ing through a series of equilibria disrupted by tearing- 
like instabilities. Similarly to the case with disordered 
vorti c al bo undary forcing velocities (Rappazzo e t al] 
120071 l2008f ) in the fully nonlinear stage the system is 
highly dynamical and chaotic (and increasingly so at 
higher Reynolds numbers). For this we do not ob- 
serve secondary tearing of the current sheets as in 2D 
high-resolution simulatio ns of decaying MHD turbu- 
lence (Biskamp & Welterl 11989 ]). as now at the small 
scales fast turbulent dynamics take place. 

2. The dynamics of the Parker model do not depend 
strongly on the pattern of the velocity forcing that 
mimics photospheric motions, as far as they are con- 
stant in time (we defer the investigation of time- 
dependent boundary forcing to a future work). The 
shear forcing [eq. ([1])] applied only at the top plate is a 
very simple and ordered one-dimensional forcing. We 
have shown that the resulting dynamics are very sim- 
ilar to those developed when a more complex and dis- 
ordered "vortex-type" forcing velocity is applied. We 
conclude that the turbulent properties of the system 
are not induced by the complexity of the path that the 
footpoints follow. It is rather the system itself to be 
intrinsically turbulent, and turbulence develops as we 
continuously inject energy at the scale of photospheric 
motions 1,000 /cm). 

3. The system reaches a statistically steady state where, 
although the footpoints of the field-lines are continu- 
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ously dragged by the forcing shear, this does not in- 
duce a sheared magnetic field in the computational 
box. In fact the topology of the magnetic field is not 
a mapping of the forcing velocity field. Nonlinear in- 
teractions are able to redistribute the energy that is 
injected only at the wavenumber = 4ex also to 
perpendicular wavenumbers and to smaller wavenum- 
ber s through an MHD turbulent cascade. In physical 
space this corresponds to the magnetic field being or- 
ganized in magnetic islands, to small-scales formation 
(current sheets elongated along the axial direction) 
and to magnetic reconnection taking place. 

4. Kinetic and Magnetic energies develop an inert ial 
range where spectra exhibit a power-law behavior. 
Fluctuating magnetic energy dominates over kinetic 
energy. Spectra and integrated quantities, like ener- 
gies and total dissipative rates, have values similar to 
those obtained with a vortex-type forcing velocity. In 
particular the total dissipation rate of the same system 
simulated with different Reynolds number appear to 
overlap each other beyond Re = 200, suggesting that 
this is independent of the Reynolds number beyond a 
threshold. 

As shown in Figures [TJ [2] and [3] initially the system un- 
til time t ^ 79rA follows the linear curves ([12]) and (fT7|) . 
Up to this point the shear velocity at the top boundary 
induces a sheared magnetic field in the volume. As dis- 
cussed in § 15.11 we have introduced a perturbation mim- 
icking those naturally present in the corona. With no 
perturbation the system would relax over the resistive 
diffusive timescale tr (~ 25 ta for run A) in a saturated 
diffusive equilibrium as described in ([12]) and ([17]) . While 
the simulation presented here used a very small ampli- 
tude for the perturbation (e = 10~^^), we have performed 
shortest simulations with different values for the ampli- 
tude. As expected for higher values of e the instability 
develops sooner and for smaller values later, always fol- 
lowing the linear curves until the instability transitions to 
the nonlinea r stage. Th e more complete and systematic 
analysis of lR.omeou et al.l (|2004 l2009f ) in 2D confirms 
t his behavior. 

iDahlburg et al.l ()2009f ) have performed a similar sim- 
ulation with a lower resolution and with a fixed value 
for the perturbation and for a time interval that covers 
only the initial stage of our simulations. They in fact 
stop right after the first big dissipative peak, that in our 
Figures m [2] and [3] corresponds at t ~ 100 r^. 

Continuing the simulation, and using a higher nu- 
merical resolution, the system reaches a statistically 
steady state where magnetic energy consistently fluctu- 
ates around a mean value and the shear is not recreated 
in the topology of the orthogonal magnetic field. Their 
analysis is then limited to a transient event taking place 
only during the early stages of the dynamics, and that 
afterward does not repeat. 

As shown in iRappazzo et al.l ()2008f ) during the linear 
stage the system is able to accumulate energy well be- 
yond the average value maintained in the nonlinear stage 
only if the boundary forcing velocity satisfies the condi- 
tion that its vorticity is constant along the streamlines. 
The sheared profiles used in this paper satisfy th i s con - 
dition as well the profile used by 'Dahlbur g et al.l ()2009[ ) 
(a linear combination of 6 sheared profiles). 



These profiles are a very small subset of all the possi- 
ble forcing profiles, and while they are very useful to get 
insight into the origin of turbulence in coronal loops they 
are not representative of the disordered photospheric mo- 
tions, for which the strong stress buildup required for sec- 
ondary instability to develop does not take place. The 
significance of their conclusions is then strongly dimin- 
ished. 

Furthermore the Parker angle for this system cannot be 
defined as the relative angle between magnetic field-lines 
at which the system becomes unstable. This is not a well- 
posed definition. In fact for given initial conditions the 
angle or equivalently the time (as the linear equation ([T2]) 
and ([TT]) imply) at which the instability develops depend 
on the value of the amplitude of the perturbation that 
we add to the system. Depending on the value of the 
perturbation the Parker angle so defined is not unique. 

On the other hand in the fully nonlinear stage the av- 
erage magnetic field line magnitude fiuctuates around a 
mean value. It is then possible to give a unique value for 
the Parker angle, defined now as the average inclination 
of the magnetic field-lines respect to the axial direction 
as done in Rappazzo et al. (2007, 2008), and as originally 
introduced by Parker (1988). 

As sum marized in ^ [4] the on e -point closure model de- 
veloped by iHevvaerts fc PriestI ()1992f ) splits the domain 
into large and small scales. They conjecture that the 
large-scale fields evolve into a stationary laminar regime, 
the field magnitudes determined by the effective diffu- 
sion coefficients. These laminar regimes correspond to 
our linear saturated diffusive regimes computed in § O 
In Figure [8] the dotted lines show such diffusive curves 
for different values of the Reynolds numbers. 

In their model the large-scale fields computed in this 
way are used to obtain the energy fiowing into the 
system for unit time at the boundary (the power) due 
do the work done by photospheric motions on the mag- 
netic field-lines footpoints. They also calculate, through 
an EDQNM approximation, the value of the spectral en- 
ergy fiux e fiowing along the inertial range at the small 
scales. Both S and e are functions of the effective diffu- 
sion coefficients, and the solution of the problem results 
requiring balance between the two powers S = e [S and e 
have both the dimension of a power, i.e. energy over time, 
as S is the Poynting flux integrated over the boundary 
surfac e and e is integrated over the whole volume as in 
Rappaz zo et all (^07)]. 

As shown in our simulations the large-scale flelds are 
not laminar, and they are st ationary only statis t ically . 
Nevertheless it is useful to use IHevvaerts fc PriestI (jl992r ) 
model in order to understand why it is not applica- 
ble. From Figure [8] we can estimate that the effective 
Reynolds number for which the diffusive regime dissipa- 
tion matches the dissipation of the simulated system is 
Reff = 150. Unfortunately this values is too low for 
their model to work. In fact for R = 150 the dynamics 
are so diffusive that only a few modes of the order of 
the injection scale (~ 1,000 /cm) are not suppressed but 
only reduced in magnitude. Therefore there is no flux of 
energy at the small scales 6 = 0. 

At a more fundamental level the idea to split the do- 
main into large and small scales does not work because 
nonlinearity cannot be conflned only at the small-scales. 
As shown by our simulations nonlinearity is at work at all 
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scales, and unfortunately this fundamental aspect cannot 
be circumvented. 

Finally the use of RMHD equations is valid as far as 
the magnetic field fiuctuations are small compared 
to the axial magnetic field Bq. This seems particularly 
apt to describe the dynamics of long-lived slender loops 
that apparently show no dynamics while shining bright 
at the resolution scale (~ 800 km) of current state-of-the- 
art X-ray and EUV imagers onboard Hinode and Stereo. 
Clearly these results do not apply to highly dynamical 
active regions where dynamics cannot be modeled as fluc- 
tuations about an equilibrium configuration. 

The series of simulations that we have performed 
proves that dragging the footpoints of magnetic field- 
lines in the Parker problem quickly triggers nonlinear 
dynamics for small values of the orthogonal magnetic 
fields, and that these small magnetic field fiuctuations 
are able to transport a considerable amount of energy to- 
ward the small scal es with the overah ener gy fiux ^ 1.6 x 
10^ erg cm~'^ CRappazzo et al.l 1200 8^ in the lower 
range of the observed constraint ^ 10^ erg cm~^ s~^. 

This prevents the orthogonal magnetic fiuctuations to 
grow to an arbitrarily high value, self-consistently lim- 
iting the dynamics of the Parker problem to small fiuc- 
tuations if the initial conditions are given by a uniform 
strong axial magnetic field. 
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